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—  >  We  discuss  shock-generated  turbulence  as  related  to  the  dynamics  of  hot  gaseous  channels 
produced  by  las^r  pulses  and  electric  discharges  in  gaseous  atmospheres.  Accurate  models  of 
channel  development  must  account  for  compressible  flows,  which  are  associated  with  shocks  that 
are  present  at  early  times,  as  well  as  the  incompressible  residual  motion  which  is  responsible  for 
channel  cooling.  We  present  an  efficient  numerical  technique  which  uses  the  flux-corrected 
transport  algorithm  and  permits  automatic  adjustment  of  the  time  step  to  account  for  both  types 
of  behavior.  Comparisons  of  two-dimensional  calculations  to  theory  and  to  experimental  data  — 
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on  laser  channels  show  good  agreement.  Our  numerical  results  reveal  a  mechanism  by  which  the 
distribution  of  turbulent  scale  lengths  is  generated  during  energy  deposition. 
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MULTIDIMENSIONAL  NUMERICAL  METHODS  FOR  THE  SIMULATION 


OF  SHOCK-GENERATED  TURBULENCE 

I.  THE  PHYSICS  OF  HOT  GASEOUS  CHANNELS 

Recent  experimental  studies  of  electric  discharges  in  air  and  laser  pulses  in  air 
and  pure  nitrogen  [1-2]  have  revealed  the  detailed  dynamics  of  the  hot  gaseous 
channels  which  remain.  The  production  and  cooling  of  such  hot,  gaseous  channels  is 
of  interest  not  only  in  laser  and  discharge  experiments,  hut  also  in  research  on  the 
physics  of  lightning  and  electron  beams.  During  and  immediately  af*er  energy 
deposition,  the  resulting  hot  gas  quickly  expands  to  achieve  pressure  equilibrium 
with  the  surrounding  atmosphere,  producing  a  shock  wave  which  propagates  away  in  a 
short  time.  The  channel  then  cools  on  time  scales  which  are  orders  of  nagnitude 
faster  than  those  characterizing  classical  (nonturbulent)  thermal  conduction-.  As  a 
channel  cools,  the  radius  increases  according  to  the  equation 

R2(t)  «  R2(t)  +  U«  (t-t),  (1) 

where  R  is  the  radius  of  the  channel,  t  is  time  measured  from  the  beginning  of  a  dis¬ 
charge  or  pulse,  t*r  is  the  time  at  which  pressure  equilibrium  is  reached,  and  a  is 
the  thermal  diffusivity.  For  electric  discharges  which  deposit  •  300-600  J/m,  measure¬ 
ments  at  early  times  give  a  •  500  cm2/s  [l].  The  experimental  C02  laser  pulses  in 
nitrogen  deposit  ■  9  J/a,  giving  a  ■  250  cm2/s  [2].  For  nonturbulent  thermal 
conduction,  we  have  a  *  1.0  cm2/s  for  air  at  300E  and  1  atm  [3 ] . 

As  a  channel  cools,  turbulent  structure  appears  first  at  the  boundary  and  then  in 
the  interior.  This  has  led  us  naturally  to  search  for  mechanisms  in  which  long-lived 
rotational  motion  is  produced.  The  only  successful  analysis  to  date  [  1 , L. ]  relies  on 
the  (rigorous)  equation  for  the  time  development  of  vorticity,  i.e.. 
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St  +  £  *  X  -  5  *fic  +  <2P  *  F>/P2  .  (2) 

where  |  *  £  x  X  is  'the  vorticity,  v  is  the  fluid  velocity,  p  is  the  density, 
and  P  is  the  pressure.  According  to  eq.  (2),  any  misalignment  between  the  pressure 
and  density  gradients  during  expansion  of  the  gas  to  pressure  equilibrium  will 
generate  vorticity.  This  will  occur  when  asymmetries  exist  in  the  structure  of  a 
pulse  and  when  consecutive  pulses  sure  noncollinear.  After  pressure  equilibrium  is 
reached,  we  may  represent  the  flow  field  in  terms  of  one  or  more  pairs  of  vortex 
filaments  of  strength  ±k,  where 

k  ■  0ffl  (Rl-  RQ)  In  (pj'po)  f.  (3) 

In  eq.  (3),  Ug  is  a  characteristic  expansion  velocity;  RQ  is  the  initial  channel 
radius;  R|  is  the  radius  of  the  channel  Just  after  expansion  to  pressure  equilibrium; 
P«  is  the  ambient  density;  p0  is  the  density  at  the  center  of  the  channel;  and 
f(<  1)  is  a  form  factor.  The  number  of  vortex  filament  pairs  depends  on  the 
situation. 

II.  NUMERICAL  MODEL 

An  accurate  numerical  model  of  channel  dynamics  is  necessary  for  the  solution  of 
problems  which  are  too  difficult  to  treat  analytically,  as  well  as  for  the  verifica¬ 
tion  and  calibration  of  our  theoretical  results  [1,1+1 .  We  favor  the  use  of  a  finite 
difference  scheme  to  Integrate  the  equations  for  conservation  of  mass,  momentum,  and 
energjr.  A  major  complication  is  the  fact  that  the  flows  vhich  generate  vorticity  are 
supersonic  and  require  short  time  steps,  as  dictated  by  the  Courant  condition,  while 
the  (incompressible)  residual  flows,  vhich  are  responsible  for  channel  cooling,  could 
be  computed  with  much  longer  time  steps  (<  100  x).  The  model  vhich  we  describe  below 
includes  a  simple,  effective  method  of  numerically  separating  the  rotational  flows 
from  the  (radial)  shock  flews  at  times  when  the  latter  are  no  longer  important  near 


the  channel.  Ibis  permits  the  adjustment  of  the  time  step  to  account  for  the  3 lover 
flews,  saving  considerable  conputer  time. 

Our  calculations  are  multidimensional,  using  time-step  splitting  in  conjunction 
with  the  latest  version  of  flux-corrected  transport  (FCT)  1 5)  •  The  grid  is  Cartesian 
and  is  uniform  and  finely  spaced  in  the  region  near  the  channel.  Outside  of  this 
fine  grid,  the  cells  increase  geometrically  so  that  the  boundaries  are  far  from  the 
channel.  This  permits  the  shock  to  decouple  completely  from  the  channel  before 
interacting  with  the  boundary.  After  the  shock  has  moved  veil  away  (>  10  channel 
radii)  from  the  channel,  the  pressure  field  near  and  vlthin  the  channel  is 
approximately  ambient.  At  this  point,  ve  calculate  and  subtract  the  average  radial 
velocity  field  from  the  total  velocity  field,  leaving  only  the  rotational  flovs  of 
interest.  To  conpute  this  average  radial  velocity  field,  ve  use  the  following  self- 
consistent  method,  based  on  digital  miltichannel  analysis: 

(1)  Relative  to  the  center  of  the  channel,  define  a  set  of  contiguous  radial 
intervals  which  enconpasa  the  entire  differencing  grid,  as  shown  in  fig.  1. 

(2)  For  each  interval,  determine  the  grid  points  whose  radial  displacements  from  the 
channel  center  fall  within  the  interval.  Average  the  radial  velocity  conponents  of 
these  grid  points. 

(3)  From  thi3  average  radial  velocity  field,  denoted  <vr>,  determine  the  average 

radial  potential  function  by  integrating  the  Poisson  equation, 

2 

7  4  s  7  <v  >. 
r  r 

(4)  Interpolate  the  field  <vr>  onto  the  grid  by  straightforward 
differencing  of  the  equation  7r  $  *  <vr>  to  maintain  consistency  with  the 
differencing  solution  for  the  total  velocity  field. 

(5)  Subtract  the  field  <vr>  from  the  total  velocity  field  at  each  grid  point. 

(6)  Set  the  pressure  field  to  the  ambient  value  and  set  the  velocities  near  and  at 
the  shock  wave  to  zero  to  eliminate  all  high  speed  flovs  and  to  prevent  shock 
interaction  with  the  boundary. 


Immediately  after  the  above  sequence,  ve  suet  increase  the  size  of  the  time  step, 
since  the  flovs  of  interest  are  ouch  slower  than  the  sound  speed.  Our  use  of  FCT, 
however,  requires  that  the  Courant  condition  be  satisfied.  Ve  accomplish  this  by 
scaling  the  pressure  down  by  a  factor  0,  so  that  P  ♦  P/8  and  cg  ♦  cg/8l/2»  where 
the  local  speed  of  sound  is  cg  ■  (yP/p)1/2  and  y  is  the  ratio  of  principal 
specific  heats,  lhe  time  step  will  scale  roughly  as  At  ♦  81/2  At.  O'Rourke  and 
Bracco  [6]  have  employed  a  similar  scaling  technique  to  reduce  the  computational  time 
required  to  model  unsteady  laminar  flames  using  multidimensional,  fully  compressible 
computer  codes.  Da  their  method,  the  speed  of  sound  remains  the  same  while  the 
velocity  and  other  related  parameters  are  scaled  to  larger  values.  Hie  transfor¬ 
mation  of  O'Rourke  and  Bracco  thus  artifically  increases  the  speed  of  the  flow 
relative  to  the  speed  of  sound  while  permitting  a  time  step  which  is  approximately 
the  same  as  that  which  would  have  been  dictated  by  the  Courant  condition  prior  to 
scaling.  Both  the  scaling  method  used  by  us  and  that  of  O'Rourke  and  Bracco 
increase  the  Mach  number  of  the  flow  to  reduce  the  effect  of  the  Courant  condition 
on  total  confuting  time. 

An  appropriate  value  of  8  will  insure  that  the  errors  introduced  by  our 

technique  are  small.  Setting  the  pressure  field  to  the  ambient  value  in  step  (6) 

above  will  initially  eliminate  the  minute  pressure  gradients  which  maintain  the 

rotational  motion  of  the  fluid.  This  will  permit  the  velocity  field  in  the  vicinity 

of  a  vortex  filament  to  expand  momentarily.  The  expansion,  however,  reduces  the 

local  density  and  pressure  approximately  adiabatically  until  the  pressure  gradients 

again  counteract  the  centrifugal  force  experienced  by  a  given  fluid  element.  In  the 

Appendix  we  show  that  the  effect  on  the  vorticity  5  is  of  order  v2/c'g2,  where 
c'  s  c_/31/2  is  the  speed  of  sound  after  scaling.  In  practice  we  choose  8  so  that 

v/c"<  1/5  to  insure  that  such  errors  are  3mall.  Another  important  concern  is  the 
effect  of  scaling  on  the  incompressibility  of  the  residual  flow  field.  Since  errors 
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related  to  compressible  effects  for  a  subsonic  flow  field  are  also  of  order  v2/c'82 
17],  our  constraint  on  0  again  insures  that  such  errors  are  less  than  a  few  percent. 
Should  0  be  too  large,  we  would  expect  some  energy  to  be  diverted  to  the  formation  of 
compressional  wa^es,  resulting  in  a  reduction  in  intensity  of  the  rotational  motion. 
Our  calculations  have  verified  that  the  flow  pattern  remains  rotationally  stable  and 
accurate  for  several  siaulated  milliseconds  after  our  procedure  is  performed.  The 
residual  vorticity  has  proven  to  be  correct  to  first  order,  as  expected,  and  no 
congressional  waves  have  arisen  at  late  times,  again  supporting  our  analysis. 

III.  NUMERICAL  SIMULATIONS 

All  numerical  simulations  to  date  have  used  a  Cartesian  grid  of  100  x  100  cells. 
Hie  central  region,  where  energy  i3  deposited,  consists  of  50  x  50  square  cells  chosen 
so  that  an  individual  pulse  would  be  contained  in  an  area  no  smaller  than  8x8 
cells.  Outside  the  uniform  central  region,  each  cell  dimension  increases 
geometrically  by  a  factor  (l+d)  from  cell  to  cell,  where  5<0.3,  so  that  the  grid 
boundary  is  far  from  the  channel.  All  pulses  initially  have  a  Bennett  pressure 
profile 

k-Eoi2  2 

P(r)  »  P„  +  (P0-P.)/(l+  — p - )  ,  W 

where  PQ  is  the  pressure  at  the  center  of  the  pulse,  P„  is  the  ambient  pressure, 
r0  is  the  position  of  the  center  of  the  pulse,  and  a  is  a  constant  (the  "Ben¬ 
nett  radius”),  -he  model  deposits  all  energy  instantaneously  in  the  form  of  an  in¬ 
crease  in  internal  energy.  We  use  a  real  air  equation  of  state  routine  based  on  the 
data  of  Gilmore  [8-9]. 

As  our  first  example,  we  choose  the  simple  case  of  noncollinear  pulses  produced 
at  different  times.  The  first  pulse  deposits  energy  at  time  t=0.0  in  the  center  of 
the  grid,  and  the  second  pulse  occurs  1  ms  later  and  is  displaced  to  the  (reader's) 
right  of  the  first  by  a  distance  equal  to  the  Bennett  radius  (a*0.5cm).  Ihe  initial 
peak  overpressures  (P0-P„)  of  the  pulses  are  4.7  atm  and  2.b  atm,  respectively. 


Figure  2  shows  density  contours  and  velocity  vectors  in  the  finely  gridded  region  at 
t*1.26  ms  and  density  contours  at  t»2ms.  The  contours  range  from  3«0xlQ“4  g/um3  at 
the  center  of  the  grid  to  l.lxlO"*3g/cm3  at  the  edge.  Fbom  Fig.  2a  we  see  that  the 
noncollinearity  of  the  pulses  produces  a  noncircular  channel  cross  section.  The 
velocity  vector  plot  clearly  demonstrates  the  presence  of  a  vortex  filament  pair 
vhose  flow  pattern  pulls  fluid  from  the  position  of  the  first  pulse  toward  that  of 
the  second.  Figure  2c  shows  that  the  gas  at  the  center  of  the  system  has  cooled 

under  the  influence  of  the  residual  flows.  These  features  are  in  complete  agreement 

with  theoretical  predictions  [U]  based  on  Eq.  (2)  and  (3). 

A  more  important  and  difficult  case  is  that  of  a  laser  pulse  with  an  approxi- 
mately  circular  envelope.  Although  turbulent  cooling  might  appear  to  be  less  im¬ 
portant  in  this  situation,  experimental  data  for  a  laser  pulse  in  nitrogen  at  1.2  atm 
give  a  ■  250cm2 /s.  laboratory  burn  patterns  have  also  revealed  that  nonuniformities 
("hot  spots")  exist  in  the  interior  of  the  laser  pulse.  We  have  used  this  fact  in  a 

test  calculation  in  which  the  actual  laser  pulse  is  assumed  to  consist  of  3even 

smaller,  identical,  simultaneous  pulses,  as  in  the  pressure  contour  diagram  of  Fig. 

3a.  Again  we  show  only  the  uniform  central  region  of  the  grid,  which  is  1.2  cm  x  1.2 
cm.  The  energy  deposited  per  unit  length  is  «  9  J/m,  in  agreement  with  experiment  [2] 
We  use  the  real  air  equation  of  state  routine  (rather  than  one  for  nitrogen)  and  1.0 
atm  ambient  pressure.  The  pressure  contours  range  from  l.lxlO6  dyne/cm2  to  3. 0x10 5 
dyne/cm2  while  the  density  contours  in  the  remaining  diagrams  range  from  5.7xl0-4 
g/cm3  to  l.lxl0"3  g/cm3.  The  rather  mild  hot  spots  in  Fig.  3a  would  be  marginally 
detectable  by  a  burn  pattern.  From  Fig.  3a,  we  see  that  the  shock  wave  produced  by 

each  hot  spot  will  sweeo  through  the  density  minima  occurring  at  the  positions  of  the 
other  hot  spots,  producing  vorticity  according  to  Eq.  (2).  By  t=56ys,  the  channel 

has  achieved  pressure  equilibrium  and  a  temperature  of  <  650K.  The  envelope 

is  nearly  circular,  even  though  nonunifcrmities  persist  in  the  interior  of  the 

channel.  Eecause  the  grid  i3  rather  coarse  and  does  not  have  sixfold  symmetry  and 


because  we  use  time-step  splitting,  the  channel  properties  do  not  retain  the  exact 
sixfold  symmetry  of  Fig.  3a.  The  use  of  a  finer  grid,  and  perhaps  multidimensional 
FCT  [lOj,  would  inprove  the  situation,  although  nultidiaensional  FCT  is  slower  (*2x) 
and  requires  more  storage  space,  as  compared  to  the  algorithm  used  here.  At  later 
times,  the  boundaries  distort  under  the  influence  of  local  vortex  filament  pairs,  and 
the  channel  dimensions  increase.  When  viewed  from  a  direction  perpendicular  to  the 
channel  axis,  the  channel  should  appear  to  expand,  and  the  boundary  distortions 
should  appear  as  striations  parallel  to  the  channel  axis.  Ihese  features  are  in  fact 
present  in  Schlieren  photographs.  CUr  calculation  yields  a  «■  If  cm2/s,  which  is 
within  a  factor  of  2.5  of  the  experimental  value  [2]. 


IV.  CONCLUDING  REMARKS 

This  method  for  calculating  the  generation  of  turbulence  by  shocks  and  the 
subsequent  cooling  of  hot  gaseous  channels  has  proven  to  be  accurate  and  efficient. 
We  have  reproduced  theoretical  predictions  for  consecutive,  nor.collinear  pulses  and 
have  verified  features  observed  experimentally  for  laser  channels  in  air.  Viewing 
the  diagrams  in  Fig.  3,  we  find  that  the  scale  lengths  of  nonuniformities  within  the 
laser  pulse  will  determine  the  scale  lengths  appearing  in  the  turbulent  structure. 
■Hie  resolution  provided  by  the  numerical  grid,  however,  places  a  lower  limit  on  the 
distribution  and  could  affect  the  interactions  between  turbulent  structures  of 
various  sizes. 
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(a)  1.26  ms  (b)  1.26  ms  (c)  2.0  ms 

Fig.  2  —  Density  contour  diagrams  and  velocity  vector  plot  for  the  case  of  two  noncollinear,  nonsimultaneous 
pulses.  The  first  pulse  occurs  at  t  »  0.0  /tis  and  the  second  at  t  —  1.0  ms.  The  appropriate  time  appears  below 
each  diagram. 


(a)  0.0  ijls 


(b)  56  fi s 


(c)  1.24  ms 


Fig.  3  —  (a)  Pressure  contours  and  (b),  (c)  density  contours  for  the  case  of  a  laser  pulse  with  a  circular  envelope 
and  a  nonuniform  interior.  Energy  deposition  occurs  at  t  ™  0.0  /us.  The  appropriate  lime  appears  below  each 
diagram. 
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APPENDIX 

Here  we  discuss  in  detail  the  errors  incurred  by  our  procedure  for 
eliminating  the  shock  wave  in  our  simulation  and  by  our  subsequent  scaling 
of  the  pressure  to  reduce  the  speed  of  sound,  thereby  increasing  the  size 
of  the  time  step  permitted  by  the  Courant  condition.  After  the  shock  wave 
has  moved  well  away  from  the  channel,  we  say  represent  the  residual  flow 
field  near  the  channel  as  a  superposition  of  vortex  filaments.  For  our 
estimate,  we  consider  an  isolated,  cylindrically  symmetric  vortex  filament 
of  finite  radius  A.  By  definition  the  vortieity  go  for  r  >  A,  where  r  is 
the  radial  cylindrical  coordinate  and  the  z-axis  coincides  with  the  axis 
of  the  filament.  Ihe  vortex  filament  strength  (or  the  circulation)  is 


k(R)  3  /  |  •  da  (1) 
S(R) 

in  which  S(R)  is  the  cross-sectional  area  of  a  cylinder  centered  on  the  z- 
axis  and  having  a  radius  R.  From  eq.  (l)  and  the  definition  above,  ic(R)  * 
k(A)  for  R  >  A.  Ihe  velocity  field  is  azimithally  symmetric. 


v(R)  «r<R)e^  -|^S^  * 


Before  our  techniques  are  iaplemented,  the  velocity  field  of  the  vortex  is 
directly  related  to  the  pressure  through  the  radial  acceleration  ar 
experienced  by  a  fluid  element  at  radius  R, 

r  *  -  s  * »  V-  <3! 

where  we  have  suppressed  the  functional  dependence  of  the  variables  on  R. 
For  the  cases  of  interest,  we  have  pv2  «  P,  so  that  the  speed  of  a  fluid 
element  is  small  compared  to  the  sound  speed  ca«  For  simplicity,  we 
could  assume  that  £  is  constant  and  directed  along  the  z-axls.  Ihen  Eq. 


(l)  becomes 


*(R)  -  itR2C 


and  Eq.  (2)  gives  us 


R  <  A 


v-  < 


5A2 
2R  ’ 


R  >  A 


(5) 


Bov  ve  employ  our  procedure  for  eliminating  the  shock  wave  and  increasing 
the  size  of  the  time  step.  1b  do  so  ve  set  the  entire  pressure  field  to 
the  ambient  value  and  scale  the  pressure  by  a  factor  8,  so  that 


p  ♦  P.  *  p!  *  p Jt' 


(6) 


After  this  occurs,  the  mLnute  pressure  gradient  in  eq.  (3)  no  longer 
exists,  and  the  fluid  will  consequently  expand.  Expansion  vill  continue 
until  pressure  gradients  satisfying  eq.  (3)  again  exist.  Ihe  fluid 
variables  vill  have  undergone  the  following  transformation: 

?'  =  ?'  ♦  p»*  <  p» 

m  m 

p  ♦  p*  <  p 


v  ♦  v*  <  v 


(7) 


Z  *  5*  <  K 

in  the  region  where  the  vorticity  is  nonnegllgihle.  In  eq.  (7)t  the 
asterisk  (*)  signifies  the  value  after  expansion,  and  all  variables  after 
expansion  depend  on  R.  The  density  and  pressure  vill  quickly  approach 
ambient  values  outside  of  the  radius  A* ( >A )  at  vhich  the  vorticity  £* 


2 


falls  to  a  negligible  value.  For  our  rather  crude  estimate,  we  average 
the  variables  over  the  region  Re  (0,A*)  to  obtain  values  representative  of 
the  interior  of  the  expanded  vortex  filament.  In  the  paragraph  below,  any 
variable  merited  by  an  asterisk  will  denote  this  average  value. 

The  proper  choice  of  the  scaling  variable  8  will  insure  that 
PJ.  »  P*  v*2,  and  the  expansion  will  be  brief  and  approximately  adiabatic. 
We  then  have 


P'»  -  p*  -  <5P' 

OB 

p*  -  p  -  5p 

(8) 

v*  ~  v  -  5v 

«•  *  5  -  55  , 

where  the  quantities  SP*",  fip,  6v,  are  all  smll  quantities.  To  confute 
the  effects  of  the  expansion,  we  use  the  adiabatic  gas  law: 


£11 

PI 


<f>T' 


(9) 


where  Y'  depends  on  and  p  and  is  approximately  constant  during  the 

expansion.  Using  eq.  (8)  in  eq.  (9),  we  obtain 

Pf  -  <5Pf 

_ (£=«£)  T'  „  x  _  lllP  . 

PI  P  ;  P 


(10) 


Using  the  expression 


Pt  _  p,. 


V  P  ~ 
r 


A* 


<5P' 

A* 


(11) 


in  eq.  (3)  with  R»A*,  where  the  centrifugal  force  is  a  maximum  (see  eq. 
(5))»  we  obtain 

p*v*2  ~  5P'  (12) 


Sotice  that  the  variables  before  scaling  also  satisfy  a  similar 
relation, 

pv2  -  5P  (13) 
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